Exact results for spatial decay of the one-body density matrix 
in low-dimensional insulators 
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We provide a tight-binding model of insulator, for which we derive an exact analytic form of the 
one-body density matrix and its large-distance asymptotics in dimensions D = 1,2. The system is 
built out of a band of single-particle orbitals in a periodic potential. Breaking of the translational 
symmetry of the system results in two bands, separated by a direct gap whose width is propor- 
tional to the unique energy parameter of the model. The form of the decay is a power law times 
an exponential. We determine the power in the power law and the correlation length in the expo- 
nential, versus the lattice direction, the direct-gap width, and the lattice dimension. In particular, 
the obtained exact formulae imply that in the diagonal direction of the square lattice the inverse 
correlation length vanishes linearly with the vanishing gap, while in non-diagonal directions, the 
linear scaling is replaced by the square root one. Independently of direction, for sufficiently large 
gaps the inverse correlation length grows logarithmically with the gap width. 
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The rapid progress in computational techniques used 
for calculating properties of solids enabled researchers to 
implement a localized real-space approach to describing 
such properties. It has made possible large-scale calcula- 
tions based on the density functional theory in par- 
ticular calculating the electronic structure of solids by 
means of 0{N) methods 0,0. In all these methods it 
is the one-body density matrix (DM) that is the central 
quantity. Its decay rate determines the degree of locality 
of all quantities relevant for physics and is decisive for 
the speed of the involved algorithms. This is why one 
observes a growing interest in localization properties of 
DM in recent years. 

However, the first result concerned with the rate of de- 
cay of DM was published by W. Kohn as early as in 1959 
Kohn demonstrated that in a one-dimensional cen- 
trosymmetric crystal, the Wannier functions, and hence 
the DM, decay exponentially with sufficiently large dis- 
tance. This result was extended to crystals of higher 
dimensionality by J. Dcs Cloizeaux 

For the reasons explained briefly above, the problem of 
the rate of decay of DM became the one of vital impor- 
tance about forty years later. The typical questions in this 
respect are concerned with the dependence of the decay 
rate on dimensionality and energy parameters (especially 
the direct-gap width) of systems under study. First, Baer 
and Gordon [£j gave general arguments that the exponen- 
tial decay, discovered by Kohn, is valid irrespectively of 
dimensionality, and moreover, the correlation length, £, 
that characterizes the exponential decay, scales with the 
direct-gap width, 5, like £ _1 ~ \/S. Then, the question 
of the rate of decay of DM was reconsidered by Ismail- 
Beigi and Arias in a general, model-independent con- 



text, for crystals of arbitrary dimensionality, that is in 
systems described by single-particle orbitals in periodic 
potentials. They found also the exponential decay law 
in insulators of arbitrary dimensionality but questioned 
the validity of the scaling ~ y/5, obtained in [||. 
They argued that instead of the square root scaling, the 
linear one, ~ 5, holds at least in the limit of van- 
ishing direct-gap width. After that, He and Vanderbilt 
[8J discovered a power-law factor, which multiplies the 
exponential factor in the asymptotic formula for DM in 
one-dimensional crystals. All these findings created a 
situation in which it was highly desirable to construct 
a D > 1 model of insulator whose DM could be inves- 
tigated exactly in the large-distance asymptotic regime, 
and the dependence of power and exponential laws of 
decay on dimensionality and direct-gap width could be 
determined unumbiguously. 

Quite recently, an attempt at constructing such a 
model was made by Taraskin et al ||. They proposed 
a simple tight-binding model, built out of two kinds of 
single-particle orbitals at each lattice site, which they 
consider to contain the basic features of an insulator. 
They succeeded in demonstrating analytically the exis- 
tence of a power-law factor and an exponential one in all 
three dimensions D = 1,2,3. However, the result was 
obtained only for very large values of one of the energy 
parameters of the model, which is in no definite relation 
with the direct-gap width [llj (unless some additional 
assumptions about the energy parameters of the model 
are made). Moreover, the model proposed by Taraskin 
et al is not a kind of crystal studied in the papers cited 
above, since it is translation invariant. 

We have succeeded in deriving analytically the exact 
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form of DM (in terms of special functions) and its large- 
distance properties for another system, which is a crys- 
tal. The system consists of one kind of single-particle 
orbitals at each lattice site, which are subjected to an 
external periodic potential. There is the unique energy 
parameter that determines the strength of the periodic 
potential. This system exhibits two bands separated by 
a direct gap whose width is proportional to the unique 
energy parameter. We have derived the power-law and 
exponential factors in the large distance asymptotics of 
DM in dimensions D — 1,2, for small and large direct- 
gap widths, and in arbitrary direction in the D = 2-case. 

Consider the model described by the following second- 
quantized Hamiltonian 



h = 22 Uiafdi + 1 2_\. ( a t a j + 



(1) 



In the above expression, i,j represent the lattice sites of 
a D-dimensional lattice, while stands for a pair of 

nearest neighbors on this lattice. The operators {a^, a^} 
create, annihilate, respectively, a spinless fermion in a 
single-particle orbital \i), belonging to an orthonormal 
basis. The value of the external potential at site j is 
Uj. Suppose that the underlying lattice is a D- dimen- 
sional simple cubic lattice, which consists of two inter- 
penetrating sublattices (the even and odd sublattices), 
such that the nearest neighbors of a site on one sublat- 
tice belong to the other one. Then, we set Uj = Ui on 
the odd sublattice and Uj = U2 on the even one. Under 
the periodic boundary conditions the Hamiltonian Q is 
block-diagonalized by the plane wave orbitals |fc) with 
the wave vector in the first Brillouin zone of the lattice. 
Specifically, shifting the energy scale to (U\ + U^)j2 and 
expressing all the energies in the units of the transfer in- 
tegral t, we obtain the upper, Aj£, and the lower, A^, 
bands of eigenenergies: 



A±=±2J( u /2) 2 + S^±2A(fc), 



(2) 



In the wave vector k is restricted to the first Bril- 
louin zone of one of the sublattices, u = {11% — U\)/2t 
is the unique energy parameter of the model, and Sk = 
\ J2j (i j) ex P j)) stands for the structure factor. 

Without any loss of generality, we can set u > 0. The 
two bands are mutually symmetric about zero and are 
separated by the gap S = 2u. For u ^ and completely 
filled lower band, the system given by Q is an insulator 
in the sense of the band theory, which we call the chess- 
board insulator. The eigenvectors corresponding to the 
eigenvalues A fe , |fc,±), respectively, are linear combina- 
tions of the vectors |fe) and \k + tt), where \ir) stands 
for the vector whose all components are equal to tt. By 
means of these eigenvectors one can calculate the zero- 



temperature, non-diagonal elements of DM: 
1 D 

(afa i+r ) = --^2S{R.(r), if<7 r = 2m+l, 



z=i 



(ata i+r ) = --(-l) CTi ft(r), if a r = 2m, (3) 



where 

TZ(r) = (2tt)- d / dfcexp(ifcr) A _1 (fe), (4) 
Jb.z. 

and SxTZ(r) = TZ{r 1 + 1, . . . , r D ) +H(r 1 - 1, ... , r D ), and 
so on for I = 2, . . . , D. In the D-dimensional integral 
is taken over the first Brillouin zone of a sublattice and 
r^O. For the discussion that follows it is convenient to 
introduce the parameter ay = X)i=i l r 'li which amounts 
to a (noneuclidean) distance between the lattice point r 
and the origin. Note that the matrix elements of DM 
depend on lZ(r) evaluated only at the points r with oy 
even. 

In the cases studied here we have found that up to a 
coefficient independent of o r , 



TZ(r) - o- r 7 exp(-cr T ./^), 



(5) 



for sufficiently large a r , with the power 7 depending only 
on the dimensionality of our insulator and £ depending 
on the direction of r and the direct-gap width. 

Specifically, in D = 1-case, the function lZ(r) at even 
points (er r = r = 2m) can be expressed by the Legendre 
function of the second kind, Q u {x): 



, . k fcos(2mk)dk (-1)™ . w 2 x , , 



where k 2 = (1 + (u/2) 2 ) -1 . For sufficiently large to, the 
above expression can be cast in the form (JSJ, where 



r 1 sln(^/l + (u/2f + u/2), 7-I/2. 



The small and large u asymptotic behaviors of £ in @ 
read: 



r 1 u s° u /2 + ..., and r 1 



In u + u 2 — 



(8) 



The inverse correlation length given by the formulae Q, 
ijHJ, and by high- accuracy numerical calculations, plot- 
ted versus u, are compared in Fig. ^ The asymptotic 
behavior (|SJ) has been also obtained numerically in one- 
dimensional crystals with periodic potentials of period 
greater than two [Io|. 

Naturally, in D — 2-case the form of lZ(r) is more 
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FIG. 1: The inverse correlation length, versus u in D — 1 
crystal: the data obtained from high-accuracy computation of 
TZ(r) given by © - circles, exact analytical result J7J - thick 
continuous curve, the asymptotic behavior of £~ , for small 
and large u, given in JH) - dotted curves. 
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T 2 (m+ 1/2) 



2tt (2/3) 2m+1 r(m + n + l)r(m-n + l) 
xF(m+ 1/2, m+ l/2;m + n + l;-l/o) 
xF(m+l/2,m+l/2;m-n + l;-l/a), (9) 

where F(a,b; c; x) stands for the Gauss hypergeometric 
function, T(a) for the Euler T-function, and the following 
notation was used 



a = 2/3(\/l + (3 2 +f3), I3 = u/A, 
n + r2 = 2m = cr r , ri — r2 = 2n. 



(10) 



To analyse TZ(r) given by ©, we first consider the di- 
agonal direction (ri = r2,n = 0). Then, TZ(r) can be 
expressed by the Legendre function of first kind 0, ^3 , 



(-1)' 



4tt 3 /3 T 2 (m + 1/2) 



F^i(l + 2/o) . (11) 



For sufficiently large a 
the form JSJ: 



ft(r) 



2m, the above 1Z(r) assumes 
(-1)" 1 exp(-2m/£) 



27T/3 



2m 



with 



r 1 =ln(v / l+^2 +/3), 7 = 1. 



(12) 



(13) 



Note that the inverse correlation length in l|13|l can be 
obtained from that in D = 1-case by a change of en- 
ergy scale, u — > u/2. Therefore, the small and large 



u asymptotic behaviour of is given by (JSJ) with the 
suitable change of the energy scale. Consequently, in 
D = 1-crystal and in the diagonal direction of D = 2- 
crystal, the inverse correlation length scales linearly with 
the direct-gap width, ~ S, for sufficiently small 5. 

Second, we choose the direction along an axis, say the 
first axis (7-2 = 0, a r = r± = 2m = 2 n). Then, the 
function TZ(r) in jjjj transforms into [ill fl~3| : 



K(r) 



(-iy 



■Q m i (l + 2a)P m _i 



1 



(14) 



with the large m asymptotics: 

m^oo (-l) m cxp(-2m/fl 

Acir) « ==■ 

tt\/2/3 2m 

a — 1\ 7r 

x cos [ m arccos | 

a + l) 4 



(15) 



where the inverse correlation length depends on u as fol- 
lows: 



r 



Tn 



(16) 



and 7 = 1. Along an axis, the first terms of the small 
and large u asymptotic expansions of £ _1 read: 



r 1 "3 

and £ _1 
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In 
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12 2 

u+- 



(17) 



Finally, we consider a general non-axial direction. Let, 
for definitness, < n < m. In this case we have not been 
able to obtain a large m asymptotics of the function TZ(r) 
in JHJ for arbitrary values of u. However, we succeeded in 
deriving the asymptotic form JJJ for sufficiently small and 
sufficiently large values of u. For sufficiently small values 
of u, we have found 7 = 1 and the inverse correlation 
length 



r X =-^ (M)m(\/T+^+/3)+ln 



_2a^" 
b + 2 V 1 b 



, (18) 



where the following notation has been used: 



b == (1 + 7i)a + ^(1 + ri) 2 a 2 + 4a?7 , 
f] = n/m= (x-l)/(x + l), X = n/r2- (19) 

Around zero, the asymptotic expansion of the above £ _1 
has the form 



r 1 u s\/^(i + li±^ 

^ V 2 ' V 24 7/ 2 



(20) 



Apparently, a straightforward way to obtain a large av 
asymptotics of TZ(r) is to expand A(fe) in powers of u~ 2 , 
for large u, carry out the integrals of products of cosine 
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FIG. 2: The inverse correlation length, £~ , versus u, for 
different directions \ in D — 2 crystal: the data obtained 
from high-accuracy computation of TZ(r) - symbols, exact 
analytical result along the diagonal 1131 and along an axis 
ill Oft - thick continuous curve, exact analytical result along a 
direction given by slope x GHt ~ thin continuous curves, exact 
analytical results for large u, given by 11711 and 112 H - dots. 
In the insert we show the log-log plot of £ _1 versus u, and the 
dotted curve represents the first term of small it asymptotics 



knowledge, we provide in this letter the first example of 
insulator being a crystal (with broken translational sym- 
metry) of dimension greater than one, where the relation 
between decay properties of DM and the direct-gap width 
can be established exactly, in the whole range of this pa- 
rameter. Specifically, for arbitrary direct-gap width the 
large-distance asymptotics of DM consists of the power 
law factor with the power 7 = D/2 and an exponential 
factor with the correlation length depending on direction 
and the direct-gap width. 

There is an excellent agreement between our analytic 
results and high-accuracy numerical results presented in 
the figures. It is worth to emphasize that, by the general 
arguments given in these results should not depend 
on the underlying lattice. 

As particularly interesting, we consider the results con- 
cerning the scaling of the correlation length £ with the 
vanishing direct-gap width 6, which is physically the most 
interesting regime, in the D = 2-crystal: ~ S in the 
diagonal direction while ~ \/~8 in non-diagonal di- 
rections. They resolve the controversy that has arisen in 
the literature in this respect 0, 0, @ : which of these two 
scalings holds in insulators? Apparently, the answer is 
more subtle than the question asked: both, depending 
on lattice direction. 

T.K. is grateful to the Institute of Theoretical Physics 
of the University of Wroclaw for kind hospitality and 
financial support. 



functions and to approximate the Euler T-functions, that 
arise, by the leading term of the Stirling's asymptotic 
expansion. Such a procedure leads to the large- distance 
asymptotics of TZ(r) (analogous to the one derived in Q), 
which consists of the power-law factor with 7 = 1 and the 
exponential factor with the inverse correlation length of 
the form 

- ]T (l + ^r'Mi+x"). (21) 

a=— 1,1 

In Fig. |2 we have displayed the analytic results for 
versus u, obtained for D = 2-crystal, and the same func- 
tion determined by means of high-accuracy numerical cal- 
culations. 

To summarize, our work has been inspired by the re- 
cent results of Taraskin et al 9], concerning the spatial 
decay of DM in a translation invariant model of insula- 
tor, where the relation between the decay and the direct- 
gap width has not been established. To the best of our 
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